Critical behavior of a water monolayer under hydrophobic confinement 
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We study by Monte Carlo simulations the low temperature phase diagram of a water monolayer 
confined between hydrophobic walls separated by h ~ 0.5 nm. By finite size scaling of the appro- 
priate order parameter, we find a liquid-liquid critical point (LLCP) in the universality class of the 
two-dimensional (2D) Ising model in the limit of infinite walls. However, for wall sizes up to hun- 
dreds of times larger than the monolayer thickness h, the LLCP is better described by the 3D Ising 
model universality class, something unexpected based on studies for simple liquids. We ascribe this 
result to the strong cooperativity and the low coordination number of the hydrogen bond network. 



PACS numbers: 64.70.Ja, 65.20.-w, 68.15.+e 

The study of nanoconfined water is of great interest 
for applications in nanotechnology and nanoscience 
The confinenient of water in quasi-one or two dimensions 
(2D) is leading to the discoveryof new and controversial 
phenomena in experiments [JQ and simulations 0, Q . 
Nanoconfinement, both in hydrophilic and hydrophobic 
materials, can keep water in the liquid phase at tem- 
peratures as low as 130 K at ambient pressure 0]. At 
these temperatures T and pressures P experiments can- 
not probe liquid water in the bulk, because water freezes 
faster then the minimum observation time of usual tech- 
niques, resulting in an experimental "no man's land" 
Nevertheless, new kind of experiments 0] and numeri- 
cal simulations can access this region, revealing interest- 
ing phenomena in the metastable state Q- In particu- 
lar, Poole et al. found, by molecular dynamics simula- 
tions of supercooled water, a liquid-liquid critical point 
(LLCP), in the "no mans land", at the end of a first- 
order liquid-liquid phase transition (LLPT) line between 
two metastable liquids phases with different density p: 
the high-p liquid at higher T and P, and the \ow-p liq- 
uid at lower T and P 0. The presence of a LLPT is 
experimentally observed in other systems Q, consistent 
with theoretical models fitted to water experimental data 
0, and is recovered by simulations of several models of 
water 0, [13, HI and other anomalous liquids 12, 13|. 
Alternative ideas, and their differences, have been dis- 
cussed and experiments on confined water in the 
"no man's land" have been seen as a way to test these 
ideas Therefore, it is relevant to understand which 
property, if any, confined water shares with bulk water. 
Several theoretical works have taken first steps in this 
direction [17|, such as, e.g., analyzing how the increasing 



of confinement affects the T and p of the water liquid-gas 
critical point [l^l . 

Here, to analyze the thermodynamic properties of wa- 
ter in confinement we consider a monolayer of water 
molecules between hydrophobic walls of area sepa- 
rated by h « 0.5 nm (Fig. [IJ. Atomistic simulations Q 
show that water under these conditions does not crys- 
tallize, but arranges in an disordered liquid layer, whose 



L 



h 



FIG. 1. Schematic view of a section of the water monolayer 
confined between hydrophobic walls of size L x L separated 
by ft ~ 0.5 nm. 



projection on one of the surfaces has square symmetry, 
with each water molecule having four nearest neighbours 
(n.n.). The molecules maximize their intcrmolecular dis- 
tance by occupying different distances from the two walls. 

We adopt a model that reproduces water properties 
[H [H [li-ill. The monolayer with i = 1, wa- 
ter molecules, each with four n.n., is partitioned into N 
equivalent cells of square section with size r = ^ LP- jN , 
equal to the average distance between water molecules, 
with V = hL'^, V/N > vq, with vq = hrl {tq = 2. 9 A) 
water excluded volume. By coarse-graining the molecules 
distance from the surfaces, we reduce our monolayer rep- 
resentation to a 2D system. We use periodic boundary 
conditions parallel to the walls to reduce finite-size ef- 
fects. We simulate constant N , P, T, allowing V{T, P) 
to change, with each cell having density pi = p{T, P) = 
N/V < pq = 1/vq. To each cell we associate a variable 
rii = irii = 1) depending if the cell i has Pi/po < 0.5 
{pi/po > 0.5). The water-water interaction is given by 



= X! - J Nub - JaN, 



coop- 



(1) 



The first term, summed over all the water molecules i and 
j at 0-0 distance , has U{r) = oo for r < tq (water 
van der Waals diameter), U{r) = 4e[(ro/r)-^^ — (rg/r)^] 
for r > Vq with e = 1.45 kJ/mol (van der Waals attraction 
energy) and U{r) = for r > rc = 25ro (cutoff). 

The second term represents the directional (cova- 
lent) component of the hydrogen bond (HB), with J = 
2.9 kJ/mol, A^HB = Z](ij) "»"j<^ff.j,ffjj number of HBs, 
with the sum over n.n., aij = 1, . . . , q bonding index of i 
to j, 6ab = 1 if a = 5, otherwise. Each water molecule 
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can form up to four HBs. A HB breaks when the OH — O 
distance > rmax ^ ^oh = 3.14A, because niUj = when 
nj > r-max = roV2 = 4.IOA (roH 0.96A). Adopting 
g = 6, we account correctly for the entropy loss due to 
HB formation and for the HB breaking if OOH > 30°, be- 
cause only 1/6 of the entire range of orientations [0, 360°] 
in the OH — O plane is associated to a bonded state. 

The third term accounts for the HB cooperativity 
due to the quantum many-body interaction |23| . with 
Jcr = 0.29 k,J/mol and iVcoop = J2i'^iJ2{i,k)i ^o-i^-xr.n 
where {I, k)i indicates each of the six different pairs of the 
four indices aij of a molecule i. To this term is due the 
0-0-0 correlation that locally leads the molecules to- 
ward an ordered configuration, that in supercooled bulk 
water is tetrahcdral at low P up to the second shell [iJ] . 
An increase of T or P partially disrupts the HB network 
and induces a more compact local structure, with smaller 
average volume per molecule. Therefore, for each HB we 
include an enthalpy increase Pwhb, where uhb/^o = 0.5 
is the average volume increase between high-p ices VI and 
VIII and low-p (tetrahcdral) ice Ih. Hence, the total vol- 
ume is y = Vb + A^HB^'HB, where Vb > Nvq is a stochas- 
tic continuous variable changing with Monte Carlo (MC) 
acceptance ru le 12 1| . Because the HBs do not affect the 
n.n. distance [2J|, we ignore their negligible effect on the 
U(r) term. Finally, we model the water-wall interaction 
by excluded volume. 

We simulate ^ 10^ state points with statistics of 
5 X 10^ independent calculations for systems with N ~ 
2500, . . . ,40000 water molecules, using a cluster MC al- 
gorithm [2lj, for a wide range of T and P, that in real 
units are as low as T = 120 K and as high as P = 0.4 GPa 
3"22|. Here, we adopt as units Ae/ks for T and 4e/i)o 
for P, with ks Boltzmann constant. 

We find (Fig. [2]) (i) a liquid-to-gas spinodal at P < 
for low T, (ii) a line of T of maximum density (TMD) 
along isobars that approaches the spinodal, without 
touching it, and continues in the locus (iii) of the line 
of T of minimum density (TminD) as in experiments 
[2^ and other models |l3l . [2q . We calculate the isother- 
mal compressibility Kt{T) = —{dh\{V) / dP)T along iso- 
bars, where (V) is the average volume, the isobaric ex- 
pansivity ap{P) = {d\n{V)/dT)p and the isobaric spe- 
cific heat Cp{P) = {d{H)/dT)p along isotherms, where 
{H) = {M') -I- P(l^) is the average enthalpy. We ver- 
ify the thermodynamic equilibrium by checking that the 
calculation of each quantity by its definition and by the 
fiuctuation-dissipation relation converge. For each quan- 
tity we find two maxima at low P. with higher-T maxima 
broader and weaker than lower-T maxima (22j . 

We find that both maxima of \ctp\ along isotherms [lo- 
cus (iv) of weak maxima |ap|^^ and (v) of stronger max- 
ima |q;p|'^] coincides within the error bars with the two 
maxima of Kt along isobars [locus (vi) of weak maxima 
and (vii) of stronger maxima if^?], consistent with 
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FIG. 2. Phase diagram of TV = 10'' water molecules in a mono- 
layer nanoconfined between hydrophobic walls. The loci (see 
text) are marked by lines (full for maxima, dashed for min- 
ima) or symbols with labels nearby or in the legend. The 
liquid-to-gas spinodal (gray line) delimits the region of sta- 
bility of the liquid with respect to the gas phase. The LLCP 
(large circle with A) is at the end of the LLPT (thick black) 
line. Error bars are of the order of the swaying of the lines or 
the size of the symbols. 



thermodynamic relations [31 . This consistency holds 
also for the locus (iv) intersecting the TMD line in its 



turning point [lj|, and for the locus (viii) of weak Cp 
maxima (Cp^) along isotherms bending toward the turn- 
ing point of the TminD line [2^ . The locus (ix) of strong 
Cp maxima (C^) coincides with loci (v) and (vii). Fi- 
nally, the loci (iv) and (vi) continue in the loci (x) |q;_p|™ 
of minima of |q!p| and (xi) if™ of minima of Kt where 
dPj dT along these loci is infinite or zero 13 , [2^ , respec- 
tively (not shown in FigH]). 

All the loci of maxima of response functions (iv)- 
(ix) converge toward the point A. Moreover, all the 
maxima along these loci increase in their values by ap- 
proaching the point A. Because the increase of response 
functions is related to the increase of fiuctuations and 
this is, in turn, related to the increase of correlation 
length ^, we calculate the spatial correlation function 
G{ra) = {(Jij{n)<Jik{n)) - {(Jij)'^ to estimate We find 
that for P below the point A, G{r) decays as an expo- 
nential with a characteristic correlation length ^. For P 
approaching the point A, G{r) is better approximated 
by a power-law decay with an exponential prefactor from 
which ^ can be extracted. At the point A, the exponential 
prefactor approaches a constant leaving the power-law as 
the dominant contribution for the decay, corresponding 
to ^ becoming of the order of the system size. We ob- 
serve that ^ has a maximum along isobars and that the 
Widom line, i.e. the locus (xii) of maxima coincides 
with the loci (v), (vii) and (ix). All these observations 
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FIG. 3. The size-dependent probability distribution Pjv for 
the rescaled order parameter x, calculated for Tc(N), Pc{N) 
and B{N) in Fig. lU has a symmetric shape that approaches 
continuously (from N = 2500, symbols at the top at a; = 0, 
to A'^ = 40000, symbols at the bottom) the limiting form 
for the 2D Ising universality class (full line), maximizing the 
difference with the 3D Ising universality class (dashed line). 
Lines connecting the symbols are guides for the eye. Error 
bars are smaller than the symbols size. 



[convergence of loci (iv)-(ix) and (xii) and the increase 
of maxima along these loci on approaching A] are con- 
sistent with the identification of A with a LLCP at the 
end of a first-order LLPT line in the P-T phase diagram 
along which the density, the energy and the entropy of the 
liquid are discontinuous, as discussed in previous works 
Next, we analyze in detail the LLCP. 

According to mixed-field finite-size sealing theory , 
a density-driven fluid-fluid phase transition is described 
by an order parameter M = p + su, where p represents 
the leading term, u is the energy density (both in in- 
ternal units) and s is the field mixing parameter. At 
the critical point the probability distribution of M is 
PnIM) (X Pdix), i.e. scales as an universal function pd, 
characteristic of the Ising fixed point in d dimensions, 
of X = B{M - A'4), where B = c^jN^^'^", {3 is the M 
critical exponent, is the ^ critical exponent, both de- 
fined by the universality class, and a a/ is a non- universal 
system-dependent parameter. We adjust B and Mc so 
that Pn^M) has zero mean and unit variance. 

By combining a set of 3 x 10"* MC simulations for ^ 300 
state points with 0.033 <T < 0.065 and O.OI <P< 0.90 
with the multiple histogram reweighting method j28| . 
and tuning s, T and P we verify that in the vicinity 
of the state point A the calculated Pn{x) has a sym- 
metric shape with respect to a; = (Fig. [3]). We find 
s = 0.25 ± 0.03 for our range of N. The resulting critical 
parameters Tc{N), Pc{N) and the normalization factor 
B{N) follow in fare agreement the expected finite-size 
behaviors with 2D Ising critical exponents (2^ (Fig. H]) . 
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FIG. 4. The size-dependent LLCP temperature Tc{N) (a) 
and pressure Pc{N) (b) (symbols), used for best-fitting Pm, 
extrapolate to Tc — 0.0597 and P^ — 0.554, respectively, fol- 
lowing the expected linear behaviors (lines), (c) The normal- 
ization factor B{N) (symbols) follows the power law function 
(dashed fine) a iV/3/<*-^. We use the d = 2 Ising critical expo- 
nents: e = 0, u and /3 = 1/8. 



From the finite-size analysis we extract the asymptotic 
values Tc = 0.0597 ± 0.0001 and Pc = 0.554 ±0.003, con- 
sistent with the state point A. However, these fits adjust 
well to the data only for large N. We, therefore, perform 
a more systematic analysis. 

For each TV, we quantify the deviation of the calculated 
p{N) from the expected p2 for the 2D Ising. Furthermore, 
due to the behavior of data for small N in Fig . |31 we cal- 
culate the deviation from the 3D Ising p^ 27[ . Following 
Liu et al. [l^l, we define the deviation from each pd as 



n Pd,peak — Pd,x=0 

with n total number of points for x, Pi{N) probability 
distribution of Xi, pd,i theoretical value for Xi, pd^peak — 
Pd,x=o difference between the distribution peak and its 
value at x = 0. 

We confirm s ~ 0.25 for p2 and find s = 0.10 ± 0.02 
for p3 for our range of N. For both W2 and W3 we 
find minima that become stronger for increasing N and 
are always close to Tc ~ 0.06 and Pc — 0.55, i.e. at 
approximately constant pc- We find that W2 decreases 
with increasing N, being as low as 0.046 for = 4 x lO'', 
with W2{N) ^ for A^ ^ 00 (Fig. E)- Therefore, for an 
infinite monolayer confined between hydrophobic walls 
separated by /i ~ 0.5 nm, the system has a LLCP that 
belongs to the 2D Ising universality class, consistent with 
our coarse-graining of the monolayer in 2D. 

However, by increasing the confinement, i.e. reducing 
A^ and L at constant p, W2 becomes larger than W3 . The 
difference W2 — W3 increases with decreasing A^, and for 
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FIG. 5. Deviations Wd of the calculated fi{N) from the d — 2 
(open symbols) and d — 3 (closed symbols) Ising universal 
function pd, as function of the inverse of the number of water 
molecules A'^ at constant p ~ pc. Lines are power-law fits. 
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FIG. 6. The free-energy cost to form an interface between 

d-l 

the two liquids coexisting at the LLCP scales as AG cc N d 
with d = 3 for AT < 10* and d = 2 for N > lO"*. 



A'' = 2500 is W3 ~ W2/3 ~ 0.05, a value approximately 
equal to the smallest W2 we find for a system ten times 
larger. In other words, by increasing the confinement of 
the monolayer at constant p, the LLCP has a behavior 
that approximates better the bulk [loj . with a crossover 
between 2D and 3D-behavior occurring at iV ~ 10^. 

This dimensional crossover is confirmed by the finite- 
size analysis of the Gibbs free energy cost AG/{kBTc) 
to form an interface between the two liquids in 
the vicinity of the LLCP, calculated as AG{N) = 
-kBTc{N)[hiPpf{x = 0)-lnPAr(a;MAx)], with P/v reach- 
ing a maximum at xmax- This quantity is expected to 
scale as AG oc iV~a~. We find that our data can be fit- 
ted as A^3 for small sizes and as TV 2 for large sizes with 
a crossover around N — 10^ (Fig. [6]). Considering the 
value of the estimated pc in real units (~ Ig/cm"^) f20| . 
the corresponding crossover wall-size is L ~ 75 nm. 

Our rationale for this dimensional crossover at fixed h 
is that, when h/L increases, the characteristic way the 
critical fluctuations spread over the system, i.e. the uni- 
versality class of the LLCP, resembles closely the bulk 
because the asymmetry among the three spatial dimen- 
sions is reduced. A similar result was found recently by 
Liu ct al. for the gas-liquid critical point of a Lcnnard- 
Jones (LJ) liquid confined between walls by fixing L and 
varying h [l8l |. However, in the case considered by Liu ct 
al. the crossover was expected because the number of lay- 
ers of particles was increased from one to several, making 
the system more similar to the isotropic 3D case. Here, 
instead, we consider always one single layer, changing the 
proportion h/L hy varying L instead of h. Therefore, it 
could be expected that the system belongs intrinsically 
to the 2D universality class for any L. 

Furthermore, the extrapolation of the results for the 
LJ liquid to the case of a monolayer with h/rQ ~ 1.7, 



as in our case, would predict a dimensional crossover at 
L/h<h [l3|- For our water monolayer, instead, wc find 
the crossover at L/h^ 150, i.e. two orders of magnitude 
larger than the L J case. This implies that the presence of 
a cooperative HB network and its low coordination num- 
ber, the main differences between water and a LJ fiuid, 
enhance drastically the spreading of the critical fluctua- 
tions along a network that at least for the flrst shell is 
similar to the one in 3D, lifting the effective dimension- 
ality of the conflned monolayer. 

In conclusion, we analyze the low-T phase diagram of a 
water monolayer confined between hydrophobic parallel 
walls of size L separated by /i « 0.5 nm. We identify 
many loci of the phase diagram, the Widom line, the 
LLPT and the LLCP. We show that the LLCP belongs 
to the 2D Ising universality class only ior L/h > 150. For 
smaller L at the same p, i.e. for stronger confinement, 
the LLCP is better described by 3D, bulk-like, critical 
behavior, as a consequence of the high coopcrativity and 
low coordination number of the HB network. 
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grant NMP4-SL-201 1-266737. 
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